romberg integration

#include <math.h>

#define ROMB_MAX 20
typedef double Value;

Value romberg(Fn f, Value a, Value b, Value eps)
{
    /* approximate integral of `f' between `a' and `b' subject
     * to a given `eps'ilon of error.
     *
     * use Romberg method with refinement substitution,
     * x = (3u-u^3)/2 which prevents endpoint evaluation and 
     * causes non-uniform sampling.
     *
     * max samples 2^(ROMB_MAX)+1 
     */
    int n, m, i;
    Value h, s;

    Value c[ROMB_MAX+1];
    Value* tmp;
    int lim;
    Value p;
    Value t, u;
    Value l;

    h = 2;
    c[0] = 0;
    l = 0;
    lim = 1;
    b -= a;
    for (n = 0; n < ROMB_MAX; ++n) {
        p = h/2 - 1;
        s = 0.0;
        for (i = 0; i < lim; i++) { 
            t = 1-p*p;
            u = p + t*p/2;
            u = (u*b + b)/2 + a;
            s = s + t*(*f)(u);
            p += h;
        }

        p = 4;      
        t = c[0];
        c[0] = (c[0] + h*s)/2.0;
        for(i = 0; i <= n; i++) {
            u = c[i+1];
            c[i+1] = (p*c[i] - t)/(p-1);
            t = u;
            p *= 4;
        }

        if (fabs(c[i] - l) < eps*fabs(c[i])*16) break;
        l = c[i];
        lim <<= 1;
        h /= 2.0;
    }
    return c[i]*b*3/4;
}